//#include <math.h>
#include <cutil.h>
#include <Particula.h>
#include <vector_types.h>

__constant__ float constGas = 5;
__constant__ float coeficienteViscosidade = 0.001;


__device__ float CalculaWme(float d2, float h, float h9) {

        
	float massaEspecifica = 0.0;
	
	float h2 = h*h;
			   		
        //CALCULO DA MASSA ESPECIFICA.
    if ( (d2 > 0) && (d2 < h2) ) 
		massaEspecifica = (315 / (64 * 3.14 * h9)) * __powf((h2 - d2), 3);
	else 
		massaEspecifica = 2.0f; //0.0;
		
	return massaEspecifica;
}

__global__ void FindNeighbors(float3* VBO, Particula* particulas, Particula* particulas2, unsigned int* values, unsigned int* keys,
        unsigned int* count, unsigned int* offSet, int numeroColunasGridCell, int numeroLinhasGridCell, const int numeroParticulas,
        float deltaX, float deltaY, float XminGrid, float YminGrid, float massaEspecificaInicial, float massa) {

    int thread = threadIdx.x;
    int block = blockIdx.x;
    int dim = blockDim.x;
    int indiceValue = dim * block + thread;

    if (indiceValue < numeroParticulas) {
	
        int posicaoParticula = values[indiceValue];
		
		particulas2[posicaoParticula].massaEspecifica = 0.0;
		
	    float h = (int)fmin (deltaX, deltaY);
		float h9 = __powf(h, 9);

        int jAtual = (VBO[posicaoParticula].x / deltaX);
        int iAtual = (VBO[posicaoParticula].y / deltaY);
        //int idCellAtual = (iAtual * numeroColunasGridCell) + jAtual;

        int iMin = (int)fmax(0, iAtual-1);
        int iMax = (int)fmin(numeroColunasGridCell - 1, iAtual+1);
		
        int jMin = (int)fmax(0, jAtual-1);
		int jMax = (int)fmin(numeroColunasGridCell - 1, jAtual+1);

        for (int i = iMin; i <= iMax; i++) 
            for (int j = jMin; j <= jMax; j++) {
			
                int idCell = (i * numeroColunasGridCell) + j;

                int kInicial = offSet[idCell];
                int kFinal = count[idCell];

                for (int k = 0; k < numeroParticulas; k++) 
				
					if ( k < kFinal ) {
	                    int indiceParticula = values[kInicial + k];

						float dx = (VBO[indiceParticula].x - VBO[posicaoParticula].x);
						float dy = (VBO[indiceParticula].y - VBO[posicaoParticula].y);
						float d2 = dx*dx + dy*dy; 
						
						particulas2[posicaoParticula].massaEspecifica += massa * CalculaWme(d2, h, h9);
					
//            	    if(k > count[idCell])
  //              	   break;
                		}
            }
    }
}

extern "C" void launch_FindNeighbors(unsigned int blocksPerGrid, unsigned int threadsPerBlock,
        float3* posicoes, Particula* particulas, Particula* particulas2, unsigned int* values, unsigned int* keys,
        unsigned int* count, unsigned int* offSet, int numeroColunasGridCell, int numeroLinhasGridCell, const int numeroParticulas,
        float deltaX, float deltaY, float XminGrid, float YminGrid, float massaEspecificaInicial, float massa) {
    // execute the kernel
    FindNeighbors <<<blocksPerGrid, threadsPerBlock>>>(posicoes, particulas, particulas2, values, keys, count, offSet, numeroColunasGridCell, numeroLinhasGridCell, numeroParticulas, deltaX, deltaY, XminGrid, YminGrid, massaEspecificaInicial, massa);
}


